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ABSTRACT 

We  describe  an  approximate  Riemann  solver  for  the  computation  of  hypervelocity  flows 
in  which  there  are  strong  shocks  and  viscous  interactions.  The  scheme  has  three  stages, 
the  first  of  which  computes  the  intermediate  states  assuming  isentropic  waves.  A  second 
stage,  based  on  the  strong  shock  relations,  may  then  be  invoked  if  the  pressure  jump  across 
either  wave  is  large.  The  third  stage  interpolates  the  interface  state  from  the  two  initial 
states  and  the  intermediate  states.  The  solver  is  used  as  part  of  a  finite-volume  code  and  is 
demonstrated  on  two  test  cases.  The  first  is  a  high  Mach  number  flow  over  a  sphere  while 
the  second  is  a  flow  over  a  slender  cone  with  an  adiabatic  boundary  layer.  In  both  cases  the 
solver  performs  well. 
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Nomenclature,  Units 

a  :  local  speed  of  sound,  m/s 

E  :  total  energy  (internal  +  kinetic),  J/kg 

e  :  specific  internal  energy,  J/kg 

h  :  specific  enthalpy,  J/kg 

M  :  Mach  number 

P  :  pressure,  Pa 

Pr  :  Prandtl  number,  (Cpft/k) 

R  :  gas  constant,  J/kg/K 

Re  :  Reynolds  number 

T  :  temperature,  I\ 

t  :  time,  s 

U  :  Riemann  invariant 

u  :  x-component  of  velocity,  m/s 

v  :  y-component  of  velocity,  m/s 

ws  :  wave  speed  used  in  the  Riemann  solver 

x  :  ^-coordinate,  m 

y  :  y-coordinate,  m 

Z  :  intermediate  variable 

a  :  weighting  function 

p  :  density,  kg/m 3 

7  :  ratio  of  specific  heats 

p  :  coefficient  of  viscosity,  Pa.s 

Superscripts 

*  :  intermediate  states  for  the  Riemann  solver 

/locally  tangent  to  the  cone  surface 

Subscripts 

MIN  :  minimum  allowable  value 

e  :  boundary-layer  edge  condition 

x ,y  :  cartesian  components 

L,  R  :  left  state,  right  state 
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1  Introduction 


In  recent  years  the  proliferation  of  relatively  fast  computers  has  popularized  the  direct  cal¬ 
culation  of  viscous,  compressible  flows  in  a  time-accurate  manner.  In  some  situations,  such 
as  the  transient  hypersonic  flow  over  a  model  in  a  shock-tunnel,  numerical  simulation  is  the 
only  way  to  extract  detailed  information  about  the  flow  field.  Such  computations  are  very 
demanding  as  there  are  both  strong  shocks  and  rarefactions  and  strong  viscous  interactions. 

This  note  describes  a  robust  Riemann  solver  for  use  in  transient  hypervelocity  flow  calcu¬ 
lations.  The  full  code  [1]  is  based  on  a  cell-centred  time-dependent  finite- volume  formulation 
of  the  axisymmetric  Navier-Stokes  equations  in  which  the  governing  equations  are  expressed 
in  integral  form  over  arbitrary  quadrilateral  cells.  The  time  rate  of  change  of  conserved 
quantities  in  each  cell  is  specified  as  a  summation  of  the  fluxes  through  the  cell  interfaces. 
The  inviscid  components  of  the  fluxes  are  computed  with  the  approximate  Riemann  solver 
while  the  viscous  fluxes  are  calculated  by  application  of  the  divergence  theorem.  At  each 
time  step,  we  first  interpolate  the  flow  state  (consisting  of  a  set  of  values  for  p,  u,  v ,  e,  P,  a) 
to  either  side  of  each  interface  at  the  start  of  the  time  step  and  then  apply  a  Riemann  solver 
to  estimate  the  flow  state  at  the  interface  during  the  time  step.  Note  that  the  solver  is 
applied  in  a  locally  rotated  frame  of  reference  in  which  the  ^-velocity  is  normal  to  the  cell 
interface. 

There  are  a  number  of  Riemann  solvers  that  can  be  used  including  “exact”  iterative 
schemes  [2]  and  approximate  (noniterative)  schemes  [3,  4].  Th-  approximate  schemes  are 
generally  less  computationally  expensive  than  the  iterative  schemes  and,  because  the  Rie¬ 
mann  solver  comsumcs  a  large  fraction  of  the  total  computational  effort,  an  approximate 
scheme  is  favoured.  Although  the  Roe-type  solver  is  popular  because  it  is  relatively  fast, 
there  are  situations  such  as  the  double-Mach-reflection  case  [5]  and  flow  over  a  sphere  [6] 
where  it  may  produce  spurious  results.  The  Osher-type  solver  [4]  is  considered  to  be  fairly 
robust  and  free  of  adjustable  parameters,  however,  we  have  experienced  some  difficulty  in 
applying  it  to  flows  with  very  strong  shocks. 

Here,  we  take  a  middle  road  between  the  fully  iterative  solvers  and  the  single-step  ap¬ 
proximate  solver  and.  at  the  cost  of  some  computational  expense,  produce  an  approximate 
solver  which  is  it  liable  in  extreme  flow  situations  and  is  vectorizable  with  current  compilers 
for  vector  computers. 
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2  Approximate  Riemann  Solver 


The  current  solver  is  a  3-stage  approximate  Riemann  solver  in  which  the  first  stage  computes 
the  intermediate  pressure  and  velocity  assuming  isentropic  wave  interaction.  A  second  stage, 
based  on  the  strong-shock  relations,  may  be  invoked  to  improve  the  first-stage  estimate  if 
the  pressure  jump  across  either  wave  are  sufficiently  large.  In  practice,  this  modification 
has  been  required  only  in  extreme  conditions  such  as  those  found  in  the  bluff-body  test  case 
(Section  3.1).  The  final  stage  is  to  select/interpolate  the  interface  state  (p,  u,  v,  e,  P,  etc) 
from  the  set  of  left ,  right  and  intermediate  states.  If  stage  2  (strong  shock  modification)  is 
not  invoked,  the  solver  is  much  like  Osher’s  approximate  Riemann  solver  [4]. 


STAGE  1:  The  first  stage  of  the  Riemann  solver  assumes  that  a  spatially  constant  left  state 
(subscript  L)  and  right  state  (subscript  R)  interact  through  a  pair  of  finite-amplitude  (and 
isentropic)  compression  or  rarefaction  waves.  Perfect  gas  relations  ([7]  cited  in  [2])  are  used 
to  obtain  the  intermediate  states  (£’,  FT)  in  the  gas  after  the  passage  of  left-moving  and 
right-moving  waves,  respectively.  The  expressions  implemented  in  the  code  are 


n  =  Pr  =  p*  =  pl 


(7-1  )(Ul-Ur) 


1  2-y/(-y  — 1) 


2ct£,(l  +  Z) 


(1) 


and 


11  j  —  —  IL 


where  the  Riemann  invariants  are 


ULZ  +  Ur 
l  +  Z 


Ul 

Ur 
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2  an 
7  ~  1 


and  the  intermediate  variable  Z  is  given  by 


OR  //^\(^-I)/(^, 

v.\Pr) 


(2) 


(3) 

(4) 


Note  that  these  expressions  involve  the  power  operator  which  is  computationally  expensive. 
For  a  limited  range  of  base  and  exponent,  the  standard  power  function  is  replaced  by  the 
approximate  expansion  [I],  In  the  exceptional  situation  o;  [U i  —  Ur)  <  0,  we  assume  that 
a  (near)  vacuum  has  formed  at  the  cell  interface  and  set  all  of  the  interface  quantities  to 
minimum  values. 


STAGE  2:  If  the  pressure  jump  across  eithe  wave  is  large  (say,  a  factor  of  10),  then  the 
guess  for  the  intei mediate  pressure  is  mod’fied  using  the  strong  shock  relations. 
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If  P*  >  10  Pi  and  P“  >  10  Pr  then  both  waves  are  taken  to  be  strong  shock  waves  and 
the  intermediate  pressure  and  .elocity  can  be  determined  directly  as 


P*  = 


1 


x  ! 


-PL 


\fPR 


and 


u  = 


L \fpR  +  sfpi 

sfpl  ul  +  \/PR  ur 


(uL  -  Ur) 


(5) 


(6) 


sJTr  +  yfPL 

If  P*  is  greater  than  Pr  or  Pr  (but  not  both),  the  stage-1  estimate  for  P*  can  be  improved 


,vith  two  Newton-Raphson  steps  of  the  form 


p*  =  P" 
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P*  <10  Pl  , 
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(9) 

-,-1 

2*r 
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P*<10Pfi  , 

P*  >  10  Pr  . 

(10) 

The  strong-shock  expressions  used  in  (9),  (10)  can  be  obtained  from  the  normal  shock  ex¬ 
pressions  in  [8]  by  taking  the  limits  as  the  pressure  jump  becomes  large.  During  the  update, 
we  ensure  that  P*  >  Pmin  where  Pmin  is  some  small  value.  After  updating  P*,  the  inter¬ 
mediate  velocity  is  evaluated  using  the  relevant  strong-shock  relation  from  (9)  or  (10). 

STAGE  3:  Now  that  we  have  computed  the  pressure  and  velocity  in  the  intermediate  regions 
behind  the  waves,  the  other  intermediate  flow  properties  may  be  evaluated.  The  interface 
conditions  used  in  the  inviscid  flux  vector  can  then  be  selected  or  interpolated  from  the  4 
flow  states  using  the  logic  shown  in  Fig.  1.  Note  that  although  only  the  left-moving  wave  is 
discussed  below,  a  similar  procedure  is  used  to  obtain  the  flow  state  behind  the  right-moving 
wave. 

If  the  pressure  rises  across  the  left-moving  wave  (i.e.  P*  >  P^,),  the  left  wave  is  assumed 
to  be  a  shock  and  density  is  obtained  from  the  Rankine-Hugoniot  relation  as 

'(7  +  l)P*  +  (7  —  1)^' 


Pl  =  PL 


[h  +  l)PL  +  (7  -\)P'\ 


(11) 


4 


The  specific  internal  energy  is  obtained  from  the  equation  of  state  as 

P* 


er  = 


(12) 


'L  (7  -  IK  ' 

and  estimates  for  the  local  speed  of  sound  (for  later  use  in  the  interpolation  of  the  interface 
properties)  are 


aL  =  \p(l  -  1  )el  ■ 

The  velocity  of  the  wave  (relative  to  the  initial  left  state)  is  given  by 


ul  ~  wsL 


7  +  7-1' 

n 


1/2 


(13) 


(14) 


2  Pl  \Pl  7  +  1 , 
where  wsi  is  the  velocity  of  the  wave  relative  to  the  cell  boundaries. 

If  the  pressure  falls  across  the  left-moving  wave  (i.e.  Pm  <  Pi ),  the  isentropic-wave 
relations  are  used  to  obtain  the  intermediate  properties.  The  local  speed  of  sound  is  obtained 
from  the  Riemann  invariant  as 


aL  =  (lrL-ui)(l  ~  l)/2  , 

while  the  specific  internal  energy  is  obtained  from  the  sound-speed  relation  as 

K)2 


e,  = 


L  (7-1)7 

The  density  is  obtained  from  the  equation  of  state  as 


P* 


Pl 


(15) 


(16) 


(17) 


(7  -  l)e! 

and  the  velocity  of  the  leading-edge  of  the  wave  (relative  to  the  initial  left  state)  is  given  by 

(18) 


ul  ~  wsl  -  aL  . 


3  Test  Cases 

3.1  High  Mach  Number  Flow  around  a  Sphere. 

The  robustness  of  the  code  is  demonstrated  by  computing  a  Mach  60  flow  over  a  7.5mm 
radius  sphere  with  a  domain  consisting  of  a  60  x  60  mesh  of  cells.  The  y  =  0  boundary. is 
the  symmetry  line  (and  stagnation  line)  while  a  tangency  condition  is  applied  at  the  surface 
of  the  sphere.  Free-stream  conditions  of 

p  =  0.5097  x  10“2  %/m3,  P  =  427.1  Pa,  e  =  2.095  x  105  J/kg, 


5 


u  =  20600  m/s,  v  =  0,  Mnomtnai  =  60.1 


are  applied  to  the  curved  inflow  boundary,  the  shape  of  which  is  derived  from  the  shock- 
position  correlations  in  [9].  blow  conditions  at  the  outflow  boundary  are  obtained  by  zero- 
order  (constant)  extrapolation.  Initial  conditions  throughout  the  domain  are  set  to 

p  =  0.5097  x  10~‘  kg/m3,  P  =  427.1  Pa,  e  =  2.095  x  105  J/kg,  u  =  0,  n  =  0. 
Despite  the  very  high  temperatures  in  the  shock  layer,  the  gas  is  considered  perfect  with 

7  =  1.4,  R  =  287  J/kg/K,  Pr  =  0.72, 

and  Sutherland’s  law  is  used  for  the  coefficient  of  viscosity.  The  Navier-Stokes  equations 
are  then  integrated  forward  in  time  using  high-order  MUSCL  interpolation  and  Euler  time¬ 
stepping  with  a  CFL  number  of  0.5. 

Figure  2  shows  the  flow  field  (pressure  and  Mach  contours)  at  t  —  13.6 ps  after  the  flow 
has  approached  steady  state.  Discrete  points  from  experimentally  derived  correlations  [9] 
are  plotted  on  the  pressure  contours.  Given  that  M  =  60  is  beyond  the  range  of  the  data 
used  for  the  correlations,  agreement  is  good.  The  largest  deviations  are  near  the  outflow 
boundary.  Profiles  of  density  and  pressure  along  the  line  of  cells  adjacent  to  the  x-axis  are 
shown  in  Fig.  3.  The  shock  appears  to  be  captured  in  2  or  3  cells  with  no  oscillation  and  the 
density  jump  is  close  to  the  ideal  strong-shock  value  of  6.  The  pressure  ratio  from  free-stream 
to  the  stagnation  [joint  is  1621  which  is  very  dose  to  the  ideal  value  of  4636  for  M  =  60  (see 
e.g.  [S]  Table  II). 

A  similar  calculation  with  M  ~  12  was  reported  in  [l]  and,  for  that  condition,  an  Osher- 
type  solver  (i.e.  stages  I  and  3  only)  failed  to  produce  a  solution.  Also,  a  finite-difference 
scheme  using  Roe-type  flux-difference  splitting  required  a  rather  large  value  for  its  entropy- 
fix  parameter  in  order  to  obtain  a  physically  reasonable  solution  (J.  White,  NASA  Langley 
Research  Centre,  private  communication).  Also  note  that,  while  the  M  =  60  shock  is  very 
strong,  the  high  temperature  in  the  region  behind  the  shock  enhances  the  viscous  dissipation 
and  may  result  in  a  smoother  solution  than  seen  at  lower  Mach  numbers. 


3.2  Flow  Over  a  Sharp  Cone 

To  illustrate  the  behaviour  of  the  solver  in  the  presence  of  strong  viscous  effects,  we  show 
the  computed  results  for  M  ~  8  flow  over  a  sharp  1°  cone.  The  axis  of  the  cone  is  aligned 
with  the  free  stream. 


6 


Two  cases  are  considered  in  which  the  cone  ilow  domain  is  discretized  as  a  set  of  100  x  60 
cells  and  100  x  90  cells.  Free-stream  conditions  of 

p  =  1 .0809  x  10“2  kg/m3,  P  =  165.51  Pa ,  e  =  3.8281  x  104  J/kg, 

T  =  53.35  K,  u  =  1164.0  m/s,  v  =  0,  Mnorninai  =  7.95 

are  applied  to  the  left  and  upper  boundaries  while  the  outflow  boundary  conditions  are 
obtained  by  extrapolation  and  the  cone  surface  is  modelled  as  a  no-slip,  adiabatic  boundary. 
To  match  the  experimental  conditions  in  [10],  the  gas  was  considered  to  be  a  perfect  gas 
with 

7  =  1.4,  R  =  287  J/kg/K,  Pr  =  0.7, 
and  viscosity  was  obtained  from  the  Sutherland  expression 

r3/2 

a  =  1.611  x  10  b  — -  Pa.s  . 

1  T+  110.33 

Based  on  free-stream  conditions  and  the  length  of  the  cone,  the  Reynolds  is  approximately 
3.3  x  10,j.  The  initial  state  of  the  flow  in  the  domain  is 

p  =  1.0809  x  10-2  kg/m3.  P  =  165.51  Pa,  e  =  3.8281  x  104  J/kg,  u  =  0,  v  =  0. 

The  Navier-Stokes  equations  are  then  integrated  forward  in  time  using  high-order  MUSCL 
interpolation  and  Euler  time-stepping  with  a  CFL  number  of  0.5. 

Figure  4  shows  the  flow  field  (pressure  and  density  contours)  t  =  22 ms  after  the  flow 
has  approached  steady  state.  The  pressure  field  is  almost  conically  symmetric,  as  per  the 
inviscid  solution  of  Taylor  and  Macoll  [11]  (see  also  [12]  Ch.  10)  and  the  shock  angle  is 
still  approximately  the  inviscid  value  of  10.5°  ([8],  Chart  4).  The  shock,  however,  is  slightly 
curved  near  the  apex  of  the  cone.  Boundary  layer  profiles  of  velocity  and  temperature  at 
x‘  —  1.0m  are  shown  in  Fig.  5  for  both  the  present  finite-volume  solutions  and  a  boundary- 
layer  solution  using  edge  conditions  of 

p,_  =  2.044  x  10~2  kg/m3,  Pe  =  416.7  Pa,  ue  =  1148.6  m/s,  Te  =  71.04  I\. 

There  is  good  agreement  between  the  present  finite-volume  solutions  and  the  spectrally- 
accurate  solution  [  1 3 j ,  especially  near  the  cone  surface.  Although  the  outer  region  of  the 
boundary  layer  is  underresolved,  (even  for  the  100  x  90  mesh)  the  finite-volume  solutions 
appear  to  be  converging  to  the  spectral  solution. 
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if  ( u *  >  0)  then 


The  contact-discontinuity  has  moved  to  the  right 
and  the  interface  state  is  determined  from  the 
L  and  L*  states. 

if  ( P *  >  PL)  then 
The  left-moving  wave  is  a  shock, 
if  (wsl  >  0)  then 

All  waves  have  moved  to  the  right. 

Interface  state  is  equal  to  L. 
else 

Interface  state  is  equal  to  L*. 
endif 
else 

The  left-moving  wave  is  a  rarefaction, 
if  (ui  —  ai  >  0)  then 
All  waves  have  moved  to  the  right. 

Interface  state  is  equal  to  L. 
elseif  -  a*L  >  0)  then 
The  rarefaction  straddles  the  interface. 
Interpolate  the  interface  state  from 
states  L  and  L* . 
else 

The  entire  rarefaction  moved  to  the 
left  of  the  interface. 

Interface  state  is  equal  to  L* . 
endif 
endif 


The  contact  discontinuity  has  moved  to  the  left 
and  the  interface  state  is  determined  from  the 
R  and  R *  states  in  a  similar  manner... 

endif 


Figure  1:  Interpolation  logic  for  the  Riemann  solver. 
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Figure  2:  M  =  60  flow  over  a  sphere  with  a  tangency  boundary  condition:  (a)  pressure 
contours;  (b)  Mach  number  contours.  denotes  the  experimental  correlation. 
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Flow  properties  for  the  cells  adjacent  to  the  (y  =  0)  stagnation  line 


y  ,  mm 


Figure  5:  Comparison  of  the  present  finite- volume  solution  with  a  spectral  solution  at  1.0m 
from  the  cone  apex:  (a)  tangential  velocity;  (b)  temperature.  Solid  line  represents  spectral 
solution,  O  denotes  100  x  60  mesh,  A  denotes  100  x  90  mesh. 
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